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Abstract 

AID  model  based  on  physical  and  electrochemical  processes  of  a  lithium  ion  cell  is  used  to  describe  constant  current  and  hybrid  pulse  power 
characterization  (HPPC)  data  from  a  6  Ah  cell  designed  for  hybrid  electric  vehicle  (HEY)  application.  An  approximate  solution  method  for  the 
diffusion  of  lithium  ions  within  active  material  particles  is  formulated  using  the  finite  element  method  and  implemented  in  the  previously  developed 
ID  electrochemical  model  as  an  explicit  difference  equation.  Reaction  current  distribution  and  redistribution  processes  occurring  during  discharge 
and  current  interrupt,  respectively,  are  driven  by  gradients  in  equilibrium  potential  that  arise  due  to  solid  diffusion  limitations.  The  model  is 
extrapolated  to  predict  voltage  response  at  discharge  rates  up  to  40  C  where  end  of  discharge  is  caused  by  negative  electrode  active  material  surface 
concentrations  near  depletion.  Simple  expressions  are  derived  from  an  analytical  solution  to  describe  solid-state  diffusion  limited  current  for  short 
duration,  high-rate  pulses. 

©  2006  Elsevier  B.V.  All  rights  reserved. 
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1.  Introduction 

Hybrid  electric  vehicles  (HE Vs)  use  a  battery  as  a  high-rate 
transient  power  source  cycled  about  a  relatively  fixed  state-of- 
charge  (SOC).  In  the  literature  however,  most  fundamentally 
based  battery  models  focus  on  predicting  energy  available  at  var¬ 
ious  constant  current  discharge  rates  beginning  from  the  fully 
charged  state  [1-3].  Cell  phone,  laptop,  and  electric  vehicle  bat¬ 
teries  are  typically  discharged  over  some  hours  and  it  is  common 
in  the  literature  to  term  discharge  rates  of  only  4  C  (four  times  the 
manufacturer’s  nominal  one  hour  Ah  rating,  lasting  on  the  order 
of  15  min)  as  “high-rate”.  In  contrast,  Hitachi  states  that  their 
5.5  Ah  HEV  cell  can  sustain  40  C  discharge  from  50%  SOC  for 
greater  than  5  s.  Phenomenological  models  capable  of  capturing 
ultra-high-rate  transient  behavior  are  needed  to  understand  and 
establish  the  operating  limitations  of  HEV  cells. 

In  the  mathematical  modeling  of  HEV  cells,  equivalent  cir¬ 
cuit  models  are  often  employed  [4-7]  and  validated  in  either 
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the  time  or  frequency  domain  [8-10].  While  the  simplicity  of 
such  models  makes  them  attractive,  unlike  fundamental  models 
they  provide  no  insight  into  underlying  physical  cell  limita¬ 
tions.  Many  good  works  do  exist  in  the  fundamental  modeling 
of  lithium  ion  cells,  though  validated  models  are  only  available 
for  cell  phone,  laptop,  and  electric  vehicle  batteries.  We  briefly 
outline  some  of  those  works. 

Doyle  et  al.  [1 1]  developed  a  ID  model  of  the  lithium  ion  cell 
using  porous  electrode  and  concentrated  solution  theories.  The 
model  is  general  enough  to  adopt  a  wide  range  of  active  mate¬ 
rials  and  electrolyte  solutions  with  variable  properties  and  has 
been  applied  in  various  studies  [1,2,12-14].  In  [1],  Doyle  et  al. 
validated  the  model  against  constant  current  data  (with  rates  up 
to  4  C)  from  similar  cells  of  three  different  electrode  thicknesses. 
Solid  and  electrolyte  phase  mass  transport  properties  were  esti¬ 
mated  to  fit  measured  data,  and  in  particular,  the  solid  diffusion 
coefficient  for  LixC6  (Ds_  =3.9  x  10“ 10  cm2  s-1)  was  chosen  to 
capture  rate  dependent  end  of  discharge.  Interfacial  resistance 
was  used  as  an  adjustable  parameter  to  improve  the  model’s  fit 
across  the  three  cell  designs.  More  recently,  the  ID  isothermal 
model  was  validated  against  a  525  mAh  Sony  cell  phone  battery 
[2].  The  authors  used  a  large  Bruggeman  exponent  correcting  for 
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Nomenclature 

as  active  surface  area  per  electrode  unit  volume 
(cm2  cm-3) 

A  electrode  plate  area  (cm2) 

A  state  matrix  in  linear  state  variable  model  state 
equation 

B  input  matrix  in  linear  state  variable  model  state 

equation 

c  concentration  of  lithium  in  a  phase  (mol  cm-3) 

C  state  matrix  in  linear  state  variable  model  output 
equation 

D  diffusion  coefficient  of  lithium  species  (cm2  s-1) 
D  input  matrix  in  linear  state  variable  model  output 
equation 

F  Faraday’s  constant  (96,487  C  mol-1) 
z'o  exchange  current  density  of  an  electrode  reaction 

(A  cm-2) 

I  applied  current  (A) 

jLl  reaction  current  resulting  in  production  or  con¬ 
sumption  of  Li  (A  cm-3) 

L  width  (cm) 

p  Bruggeman  exponent 

Q  capacity  (A  s) 

r  radial  coordinate  (cm) 

R  universal  gas  constant  (8.3143  J  mol-1  K-1) 

Rf  him  resistance  on  an  electrode  surface  (Q  cm2) 
Rs  radius  of  solid  active  material  particles  (cm) 

Rsei  solid/electrolyte  interfacial  him  resistance 

(£2  cm2) 

s  Laplace  variable  (rad  s  - 1 ) 

t  time  (s) 

t[ 3  transference  number  of  lithium  ion  with  respect 

to  the  velocity  of  solvent 
T  absolute  temperature  (K) 

Ts  time  step  (s) 

U  open-circuit  potential  of  an  electrode  reaction  (V) 

v  negative  electrode  solid  phase  stoichiometry  and 

spatial  coordinate  (cm) 

y  positive  electrode  solid  phase  stoichiometry 

Greek  symbol 

o?a,  ac  anodic  and  cathodic  transfer  coefficients  for  an 
electrode  reaction 
8  penetration  depth  (cm) 

£  volume  fraction  or  porosity  of  a  phase 

0  surface  overpotential  of  an  electrode  reaction  (V) 

k  conductivity  of  an  electrolyte  (S  cm-1) 

/cd  diffusional  conductivity  of  a  species  (A  cm-1) 

o  conductivity  of  solid  active  materials  in  an  elec¬ 

trode  (S  cm-1) 

r  dimensionless  time  for  solid-state  diffusion 

0  volume-averaged  electrical  potential  in  a  phase 

(V) 

co  frequency  (rad  s - 1 ) 


Subscripts 

e 

electrolyte  phase 

s 

solid  phase 

s,avg 

average,  or  bulk  solid  phase 

s,e 

solid  phase  at  solid/electrolyte  interface 

s,max 

solid  phase  theoretical  maximum  limit 

sep 

separator  region 

— 

negative  electrode  region 

+ 

positive  electrode  region 

Superscripts 

eff 

effective 

Li 

lithium  species 

tortuosity  in  the  negative  electrode  (p  =  33)  leading  to  the  con¬ 
clusion  that  the  battery  was  electrolyte  phase  limited.  Though 
the  model  successfully  predicted  end  of  discharge  for  rates  up 
to  3  C,  the  voltage  response  during  the  hrst  minutes  of  discharge 
did  not  match  and  was  found  to  be  sensitive  to  values  chosen  for 
interfacial  resistances. 

While  the  majority  of  the  modeling  literature  is  devoted  to 
voltage  prediction  during  quasi- steady  state  constant  current 
discharge  and  charge,  we  note  several  discussions  of  transient 
phenomena  relevant  to  HEV  cells.  Neglecting  effects  of  con¬ 
centration  dependent  properties  (that  generally  change  mod¬ 
estly  with  time),  the  three  transient  processes  occurring  in  a 
battery  are  double-layer  capacitance,  electrolyte  phase  diffu¬ 
sion,  and  solid  phase  diffusion.  Due  to  the  facile  kinetics  of 
lithium  ion  cells,  Ong  and  Newman  [15]  demonstrated  that 
double-layer  effects  occur  on  the  millisecond  time  scale  and  can 
thus  be  neglected  for  current  pulses  with  frequency  less  than 
-100  Hz. 

Unlike  double-layer  capacitance,  electrolyte  and  solid  phase 
diffusion  both  influence  low-frequency  voltage  response  and  the 
relative  importance  of  various  diffusion  coefficient  values  can 
be  judged  either  in  the  frequency  domain  [16,17]  or  through 
analysis  of  characteristic  time  scales  [13,14].  Fuller  et  al.  [14] 
studied  the  practical  consequence  of  these  transient  phenom¬ 
ena  by  modeling  the  effect  of  relaxation  periods  interspersed 
between  discharge  and  charge  cycles  of  various  lithium  ion 
cells.  Voltage  relaxation  and  the  effect  of  repeated  cycling  were 
influenced  very  little  by  electrolyte  concentration  gradients  and 
were  primarily  attributed  to  equalization  of  local  state  of  charge 
across  each  electrode.  Non-uniform  active  material  concentra¬ 
tions  would  relax  via  a  redistribution  process  driven  by  the 
corresponding  non-uniform  open-circuit  potentials  across  each 
electrode. 

This  work  extends  a  previously  developed  ID  electrochem¬ 
ical  model  [3]  to  include  transient  solid  phase  diffusion  and 
uses  it  to  describe  constant  current,  pulse  current,  and  driving 
cycle  test  data  from  a  6  Ah  lithium  ion  cell  designed  and  built 
for  the  DOE  FreedomCAR  program.  The  model  highlights  sev¬ 
eral  effects  attributable  to  solid-state  diffusion  relevant  to  pulse 
operation  of  HEV  batteries. 
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Fig.  1.  Schematic  of  ID  (x-direction)  electrochemical  cell  model  with  coupled 
ID  (r-direction)  solid  diffusion  submodel. 

2.  Model  formulation 

The  ID  lithium  ion  cell  model  depicted  in  Fig.  1  consists  of 
three  domains — the  negative  composite  electrode  (with  hixC^ 
active  material),  separator,  and  positive  composite  electrode 
(with  a  metal  oxide  active  material).  During  discharge,  lithium 
ions  inside  of  solid  LixC6  particles  diffuse  to  the  surface  where 
they  react  and  transfer  from  the  solid  phase  into  the  electrolyte 
phase.  The  positively  charged  ions  travel  via  diffusion  and 
migration  through  the  electrolyte  solution  to  the  positive  elec¬ 
trode  where  they  react  and  insert  into  metal  oxide  solid  particles. 
The  separator,  while  conductive  to  ions,  is  an  electronic  insula¬ 
tor,  thus  forcing  electrons  to  follow  an  opposite  path  through  an 
external  circuit  or  load. 

The  composite  electrodes,  consisting  of  active  material  and 
electrolyte  solution  (along  with  lesser  amounts  of  conductive 


Table  1 

Governing  equations  of  lithium  ion  cell  model 


Conservation  equations 


Charge 

Electrolyte  phase 
Solid  phase 


d 

dx 


K 


eff 


d 

dx 


0e  )  + 


d 

dx 


K 


eff 

D 


d 

dx 


In  Cq  ]  +  jLl  = 


IK* 


0  (1) 
(2) 


Species 

Electrolyte  phase 
Solid  phase 


9(£ece)  _  9 

dt  dx 
=  Ds  _9_ 
dt  r 2  dr 


dcs 

dr 


+ 


1  -t 


0 

+ 


F 


(3) 

(4) 


filler  and  binder,  not  shown  in  Fig.  1),  are  modeled  using 
porous  electrode  theory,  meaning  that  the  solid  and  electrolyte 
phases  are  treated  as  superimposed  continua  without  regard  to 
microstructure.  Electrolyte  diffusion  and  ionic  conductivity  are 
corrected  for  tortuosity  resulting  from  the  porous  structure  using 
Bruggeman  relationships,  Dff  =  De£e  and  /ceff  =  k£q,  respec¬ 
tively. 

Electronic  conductivity  is  corrected  as  a  function  of  each 
electrode’s  solid  phase  volume  fraction,  <jeff  =  <j£s.  Mathemati¬ 
cal  equations  governing  charge  and  species  conservation  in  the 
solid  and  electrolyte  phases  are  summarized  in  Table  1 . 

Distribution  of  liquid  phase  potential,  <fiQ,  is  described  by 
ionic  and  diffusional  conductivity,  Eq.  (1),  with  diffusional  con¬ 
ductivity: 


2RTKett 

y 


it 


0 


1) 


dln/±\ 

dlnce  J 


described  by  concentrated  solution  theory  [3,11].  Distribution 
of  solid  phase  potential,  0S ,  is  governed  by  Ohm’s  law,  Eq. 
(2).  The  reaction  current  density,  ju ,  sink/source  term  appears 
with  opposite  signs  in  the  charge  conservation  equations  for  the 
two  phases,  maintaining  electroneutrality  on  both  a  local  and 
global  basis.  Reaction  rate  is  coupled  to  phase  potentials  by  the 
Butler- Volmer  kinetic  expression: 


asio 


RT 


tfsEI  Li 

a& 


—exp 


RT 


^SEI  Li 

a& 


with  overpotential,  r\,  defined  as  the  difference  between  the 
solid  and  liquid  phase  potentials,  minus  the  open-circuit  poten¬ 
tial  of  the  solid  or  r\  -  0S  —  0e  —  U.  Exchange  current  den¬ 
sity,  io,  exhibits  modest  dependency  on  electrolyte  and  solid 
surface  concentrations,  ce  and  cs,e,  respectively,  according  to 
k  =  (ce)“a(cs,max  -  cSje)aa(cs,e)“c-  A  resistive  film  layer,  tfSEi, 
may  be  included  to  model  a  finite  film  at  the  surface  of  elec¬ 
trode  active  material  particles  which  reduces  the  overpotential’s 
driving  force  [18]. 

Solid-state  transport  of  Li  within  spherical  LixC6  and  metal 
oxide  active  material  particles  is  described  by  diffusion,  Eq.  (4). 
Note  that  the  macroscopic  cell  model  requires  only  the  value  of 


Boundary  conditions 


90e 

dx 


x=0 


90e 

dx 


=  0 

x=L 


.eff  90s 
—  dx 


x=0 


.eff  90s 
+  dx 


x— L_+Lsep+L-|- 


A  ’ 


90s 

dx 


x=L_ 


90s 

dx 


x=L-+L 


sep 


=  o 


dce 

dx 


=0 


dce 

dx 


=  0 


dcs 

dr 


r= 0 


=  0,  D 


deg 
s  dr 


 J 


Li 


r=RK 


asF 
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solid  phase  concentration  at  the  particle  surface,  cs?e  =  cs\r=Rs , 
to  evaluate  local  equilibrium  potential,  U,  and  exchange  current 
density,  z'o-  Authors  have  used  various  approaches  in  their  treat¬ 
ment  of  Eq.  (4),  including  an  analytical  solution  implemented  as 
a  Duhamel  superposition  integral  [1 1],  a  parabolic  concentration 
profile  model  [19],  and  a  sixth-order  polynomial  profile  model 
[20] .  In  the  present  work  we  approximate  Eq.  (4)  using  the  finite 
element  method  as  described  in  Appendix  A.  Spatial  discretiza¬ 
tion  with  five  linear  elements  unevenly  spaced  along  the  particle 
radius  provides  sufficient  resolution  of  cs,e(0  as  a  function  of 
ju(t)  at  both  short  and  long  times.  The  various  approaches  to 
solid-state  diffusion  modeling  are  contrasted  and  discussed  in 
Section  4. 

For  numerical  solution,  the  ID  macroscopic  domain  is 
discretized  into  approximately  70  control  volumes  in  the  x- 
direction.  The  solid  diffusion  submodel  (Appendix  A)  is  sep¬ 
arately  applied  within  each  control  volume  of  the  negative  and 
positive  electrodes.  The  four  governing  equations  (Table  1)  are 
solved  simultaneously  for  field  variables  ce,  cs,e,  0e,  and  0S. 
Current  is  used  as  the  model  input  and  boundary  conditions 
are  therefore  applied  galvanostatically.  Cell  terminal  voltage  is 
determined  by  the  equation: 

Rf 

V  —  0s \x=L  0s  lx=0  —I,  (7) 

A 

where  Rf  represents  a  contact  resistance  between  current  collec¬ 
tors  and  electrodes. 

3.  Model  parameterization 

Low-rate  static  discharge/charge,  hybrid  pulse  power  char¬ 
acterization  (HPPC),  and  transient  driving  cycle  data  were  pro¬ 
vided  by  the  DOE  FreedomCAR  program  for  a  276  V  nominal 
HEV  battery  pack  consisting  of  72  serially  connected  cells.  Data 
was  collected  according  to  Freedom  CAR  test  procedures  [4]. 
For  the  purpose  of  HEV  systems  integration  modeling  [21],  we 
were  tasked  to  build  a  mathematical  model  of  single  cell  of  that 
pack.  We  make  no  attempt  to  account  for  cell-to-cell  variability 
and  present  all  data  on  a  single  cell  basis  by  dividing  measured 
pack  voltage  by  72.  Due  to  the  proprietary  nature  of  the  proto¬ 
type  FreedomCAR  battery  we  were  unable  to  disassemble  cells 
to  measure  geometry,  composition,  etc.,  and  thus  adopt  values 
from  the  literature  and  adjust  them  as  necessary  to  fit  the  data. 
By  expressing  capacity  of  the  negative  and  positive  electrodes 
as 

Q—  —  £s—  A)(cs?max_)(Ax)F, 

G+  =  £s+(F+A)(cs?max+)(Ay)JF,  (8) 

low-rate  capacity  data  provides  a  rough  gauge  of  electrode 
volume  and  stoichiometry  cycling  range,  assuming  electrode 
composition  and  electrode  mass  ratio  values  from  Ref.  [2] .  This 
mass  ratio  is  later  shown  to  result  in  a  well-balanced  cell  at  both 
high  and  low  rates.  Discharge  capacity  at  the  1  C  (6  A)  rate  was 
measured  to  be  7.2  Ah  and  we  define  stoichiometry  reference 
points  for  0%  and  100%  SOC  (listed  in  Table  2)  on  a  7.2  Ah 
basis. 


Fig.  2.  Empirical  open-circuit  potential  relationships  for  negative  and  positive 
electrodes. 

The  negative  electrode  active  material  almost  certainly  con¬ 
sists  of  graphite  (LixC6)  given  its  widespread  use  in  reversible 
lithium  ion  cells.  Shown  in  Fig.  2,  we  use  the  empirical  corre¬ 
lation  for  LixC6  open-circuit  potential,  U- ,  from  Ref.  [2] .  The 
positive  electrode  active  material  could  consist  of  LiyM^CU, 
LlyCoC^,  LivNi02,  or  some  combination  of  metal  oxides.  Listed 
in  Table  2,  we  fit  our  own  correlation  for  U+  by  subtracting  U- 
from  the  cell’s  measured  open-circuit  voltage. 

Fig.  3  compares  the  model  using  parameters  listed  in  Table  2 
to  1  C  (6  A)  constant  current  discharge  and  charge  data.  This 
low-rate  data  set  is  relatively  easy  to  fit  as  it  deviates  little  from 
open-circuit  voltage. 

In  contrast,  the  voltage  perturbation  of  the  transient  HPPC 
data  set  is  more  difficult  to  fit.  The  HPPC  test  procedure,  defined 
in  [4],  consists  of  a  30  A  discharge  for  18  s,  open-circuit  relax¬ 
ation  for  32  s,  22.5  A  charge  for  10  s,  followed  by  open-circuit 
relaxation  as  shown  in  the  top  window  of  Fig.  4.  The  onset 
of  constant  current  discharge  and  charge  portions  of  the  HPPC 
profile  are  accompanied  by  brief  (0.1s)  high-rate  pulses  to 


Fig.  3.  Model  validation  versus  constant  current  charge/discharge  data  at  1  C 
(6  A)  rate. 
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Table  2 

FreedomCAR  cell  model  parameters 


Parameter 

Negative  electrode 

Separator 

Positive  electrode 

Design  specifications  (geometry  and  volume  fractions) 

Thickness,  8  (x  10-4  cm) 

50 

25.4 

36.4 

Particle  radius,  Rs  (x  10~4  cm) 

1 

1 

Active  material  volume  fraction,  es 

0.580 

0.500 

Polymer  phase  volume  fraction,  sp 

0.048 

0.5 

0.110 

Conductive  filler  volume  fraction,  £f 

0.040 

0.06 

Porosity  (electrolyte  phase  volume  fraction),  se 

0.332 

0.5 

0.330 

Solid  and  electrolyte  phase  Li+  concentration 

Maximum  solid  phase  concentration  cs?max  (x  10  3  mol  cm  3) 

16.1 

23.9 

Stoichiometry  at  0%  SOC,  xo%,  yo% 

0.126 

0.936 

Stoichiometry  at  100%  SOC,  xioo%>  yioo% 

0.676 

0.442 

Average  electrolyte  concentration,  ce  (x  10-3  mol  cm- 

3> 

1.2 

1.2 

1.2 

Kinetic  and  transport  properties 

Exchange  current  density,  i o  (x  10-3  A  cm-2) 

3.6 

2.6 

Charge-transfer  coefficients,  aa,  ac 

0.5,  0.5 

0.5,  0.5 

SEI  layer  film  resistance,  Rsei  (£2  cm2) 

0 

0 

Solid  phase  Li  diffusion  coefficient,  Ds  (x  10“ 12  cm2  s 

-1) 

2.0 

3.7 

Solid  phase  conductivity,  a  ( S  cm-1) 

1.0 

0.1 

Electrolyte  phase  Li+  diffusion  coefficient,  De  (x  10-6 

cm2  s-1) 

2.6 

2.6 

2.6 

Bruggeman  porosity  exponent,  p 

1.5 

1.5 

1.5 

Electrolyte  phase  ionic  conductivity,  k  (S  cm- 1 ) 

k  =  0.0158ce  exp(0.85Cg4) 

K  = 

0.0158ce  exp(0.85Cg4) 

tc  =  0.0158ce  exp(0.85c<!'4) 

Electrolyte  activity  coefficient,  f± 

1.0 

1.0 

1.0 

Li+  transference  number, 

0.363 

0.363 

0.363 

Parameter 

Value 

Equilibrium  potential 

Negative  electrode,  U-  (V) 

Positive  electrode,  U+  (V) 

Plate  area-specific  parameters 

Electrode  plate  area,  A  (cm2) 

Current  collector  contact  resistance,  Rf  (Q  cm2) 

U-(x)  =  8.00229  +  5.0647*  -  12.578*1/2  - 
0.46016  exp[15.0(0.06  —  x)]  —  0.55364  exp[- 
U+(y)  =  85.681/  -  357.70/  +  613.89/  - 
0.30987  exp(5.657/150)  +  13.1983 

10452 

20 
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-2.4326(v 
555. 65 y3 
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Fig.  4.  Model  validation  versus  HPPC  test  data.  SOC  labeled  on  6Ah-basis 
per  FreedomCAR  test  procedures.  SOC  initial  conditions  used  in  7.2  Ah-basis 
model  are  41.7%,  50.0%,  58.3%,  66.6%,  and  75.0%. 


estimate  high-frequency  resistance.  Unable  to  decouple  values 
of  SEI  layer  resistance  from  contact  film  resistance  (or  cell- 
to-cell  interconnect  resistance  for  that  matter),  we  fit  ohmic 
perturbation  using  a  contact  film  resistance  of  Rf  =  20  £2  cm2. 

Neglecting  double-layer  capacitance  for  reasons  noted  ear¬ 
lier,  the  only  transient  phenomena  accounted  for  in  the  gov¬ 
erning  equations  (here  and  in  other  work  [14])  are  electrolyte 
diffusion  and  solid  diffusion.  A  parametric  study  showed  that 
while  it  was  possible  to  match  the  observed  voltage  drop  at 
the  end  of  the  HPPC  30  A  discharge  by  lowering  De  several 
orders  of  magnitude  from  a  baseline  value  of  2.6  x  10  6  cm2  s  1 
[2],  the  voltage  drop  at  short  times  was  too  severe.  Signifi¬ 
cant  decrease  in  De  also  caused  predicted  voltage  to  diverge 
from  the  measured  voltage  over  time  due  to  severe  electrolyte 
concentration  gradients.  While  recent  LiPF6 -based  electrolyte 
property  measurements  [22]  show  diffusion  coefficient,  De,  and 
activity  coefficient,  f±,  both  exhibiting  moderate  concentration 
dependency,  it  is  beyond  the  scope  of  this  work  to  consider 
anything  beyond  the  first  approximation  of  constant  De  and 
unity  f±. 

In  investigating  solid-state  diffusion  transient  effects,  we 
note  that  measured  voltage  response  only  allows  observation 
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of  characteristic  time  t  =  R^/Ds  and  will  not  provide  Rs  and 
Ds  independently.  SEM  images,  such  as  that  shown  by  Dees 
et  al.  [23]  (their  Fig.  1)  of  a  LiNi0.sCo0.15Al0.05O2  compos¬ 
ite  electrode,  often  show  bulk  or  “secondary”  active  material 
particles  (with  radii  ^5  p,m)  having  finer  “primary”  particles 
(with  radii  ~0.5  jxm)  attached  to  the  surface.  Dees  achieved 
good  description  of  LiNi0.sCo0.15Al0.05O2  impedance  data  in 
the  0.01-1  Hz  frequency  range  using  a  characteristic  diffusion 
length  of  1.0  pun.  As  active  material  composition  and  structure 
are  unknown  for  the  present  cell,  we  adopt  this  value  as  the 
particle  radius  in  both  electrodes. 

Though  hixCe  is  often  reported  to  have  more  sluggish  dif¬ 
fusion  than  common  positive  electrode  active  materials,  a  para¬ 
metric  study  on  Ds_  using  the  transient  solid  diffusion  model 
found  no  value  capable  of  describing  both  the  ~0.047  V  drop 
in  cell  voltage  from  2  to  20  s  of  the  HPPC  test  as  well  as  the 
slow  voltage  relaxation  upon  open-circuit  at  20  s.  Assuming  for 
the  moment  that  the  ^0.047  V  drop  is  caused  solely  by  solid 
diffusion  limitations  in  the  negative  electrode,  we  estimate  that 
surface  concentration  cs?e-  would  need  to  fall  from  its  initial 
value  by  4.3  x  10-4  mol  cm-3  (a  substantial  amount)  to  cause 
the  observed  0.047  V  change  in  U-.  Wang  and  Srinivasan  [24] 
give  an  empirical  formula  for  the  evolution  of  a  concentration 
gradient  within  a  spherical  particle  subjected  to  constant  surface 
flux  as 


cs,e(0  cs,avg(0  —  j 


Li 


-R< 


5a*FD< 


1  —  exp  — 


20  ^/LKf 

T  R I 


(9) 


At  steady  state  (where  the  exponential  term  goes  to  zero)  and 
with  the  assumption  of  uniform  reaction  current  density,  j ^  = 
1/  AL- ,  we  manipulate  Eq.  (9)  to  obtain  a  rough  estimate  of  the 
negative  electrode  diffusion  coefficient: 


Rs  l 

5FasAcs-  AL- 


(10) 


of  1.6  x  10  12  cm2  s  1 .  A  corresponding  characteristic  time 
whereby  the  operand  of  the  exponential  term  in  Eq.  (9)  equals 
unity  is  140s.  Our  conclusion  is  that,  while  a  negative  elec¬ 
trode  solid  diffusion  coefficient  of  Ds_  =  1.6  x  10  12  cm2  s  1 
might  cause  the  cell  voltage  to  drop  ^0.047  V,  that  voltage 
drop  would  take  much  longer  to  develop  than  what  we  observe 
in  the  data.  Repeating  calculations  for  concentration  gradient 
magnitude  and  characteristic  time  under  a  variety  of  condi¬ 
tions  revealed  that  the  observed  transient  behavior  might  be 
described  by  solid-state  diffusion  in  the  negative  electrode  if 
the  slope  dU-/dcs-  were  roughly  eight  times  steeper.  This 
is  indeed  the  case  in  the  positive  electrode,  where  at  50% 
SOC  the  open-circuit  potential  function  has  almost  seven  times 
greater  slope  with  respect  to  concentration  than  the  negative 
electrode.  Chosen  via  parametric  study,  final  values  of  Ds_ 
(2.0  x  10-12  cm2  s-1)  and  Ds+  (3.7  x  10~12  cm2  s-1)  represent 
HPPC  voltage  dynamics  in  Fig.  4  quite  well,  although  we  note 
they  are  dependent  upon  our  choice  of  particle  radius.  Were  we 
to  chose  particle  radii  of  5  pun  rather  than  1  pun,  our  diffusion 
coefficient  would  be  25  times  higher  to  maintain  characteristic 


Fig.  5.  Nominal  model  compared  to  models  where  limitations  of  electrolyte 
phase,  negative  electrode  solid  phase,  and  positive  electrode  solid  phase  diffusion 
have  been  individually  (not  sequentially)  removed. 

time  t  =  R2/Ds  and  match  the  voltage  dynamics  of  the  HPPC 
test. 

Fig.  5  quantifies  voltage  polarization  resulting  from  diffu- 
sional  transport  by  individually  raising  each  diffusion  coefficient 
by  five  or  more  orders  of  magnitude  such  that  it  no  longer  affects 
cell  voltage  response.  Despite  comparable  values  of  Ds_  and 
Ds+,  the  positive  electrode  polarizes  transient  voltage  response 
more  significantly  due  to  its  stronger  open-circuit  potential  cou¬ 
pling. 

Fig.  6  compares  model  voltage  prediction  to  data  taken  on 
the  FreedomCAR  battery  whereby  an  ABC- 150  battery  tester 
was  used  to  mimic  a  power  profile  recorded  from  a  Toyota  Prius 
HEV  on  a  federal  urban  (FUDS)  driving  cycle.  Only  the  first 
150  s  are  shown,  though  results  are  representative  of  the  entire 
test. 


Fig.  6.  Model  validation  versus  transient  FUDS  cycle  HEV  data.  Power  profile 
of  data  mimics  that  recorded  from  a  Toyota  Prius  (passenger  car)  HEV  run  on  a 
chassis  dynamometer. 
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Fig.  7.  Distribution  of  reaction  current  across  cell  for  first  30  s  of  HPPC  test  at 
58.3%  SOC. 

4.  Results  and  discussion 

4.1.  Reaction  dynamics 

Despite  aforementioned  uncertainties  in  cell  design  and 
choice  of  model  parameters,  the  model  is  still  useful  in  eluci¬ 
dating  pulse  discharge  and  charge  dynamics  resulting  from  solid 
phase  transport  limitations.  As  a  basis  for  the  subsequent  discus¬ 
sion  we  use  simulation  results  from  the  first  30  s  of  the  58.3% 
SOC  HPPC  case  (whose  voltage  response  is  denoted  with  circles 
in  Fig.  4). 

Approximately  one  second  into  the  test,  a  30  A  discharge 
current  is  applied  resulting  in  the  step  change  in  local  reaction 
current,  jL\  shown  in  Fig.  7.  At  the  onset  of  the  step  change, 
solid  phase  surface  concentrations  are  uniform  and  the  initial 
distribution  of  reaction  across  each  electrode  is  governed  by  rel¬ 
ative  magnitudes  of  exchange  current  density,  electrolyte  phase 
conductivity,  and  solid  phase  conductivity.  The  solid  phase  is  a 
much  better  conductor  than  the  electrolyte  phase  and  the  reaction 
is  distributed  such  that  Li+  ions  favor  a  path  of  least  resistance, 
traveling  the  shortest  distance  possible  in  the  electrolyte  phase. 
Negative  electrode  reaction  is  less  evenly  distributed  than  pos¬ 
itive  electrode  reaction  predominantly  due  to  the  greater  solid 
phase  conductivity  of  the  negative  electrode  (cr_  =  1.0  S  cm-1 
versus  cr+  =0.1  S  cm-1). 

As  a  consequence  of  the  initial  peak  in  reaction  current  at 
the  separator  interface,  Li  surface  concentration  changes  most 
rapidly  at  that  location  in  each  electrode,  as  shown  in  Fig.  8. 
The  effect  is  more  pronounced  in  the  negative  electrode  where 
the  larger  initial  peak  in  current  density  quickly  causes  a  gra¬ 
dient  in  active  material  surface  concentration,  dcs^/dx,  to  build 
across  that  electrode.  Local  equilibrium  potential,  U- ,  falls  most 
rapidly  at  the  negative  electrode/separator  interface,  penalizing 
further  reaction  at  that  location  and  driving  the  redistribution  of 
reaction  shown  in  Fig.  7.  As  discharge  continues,  progressively 
less  reaction  occurs  at  the  separator  interface  and  more  reaction 
occurs  at  the  current  collector  interface. 

Positive  electrode  reaction  current  exhibits  similar  redistri¬ 
bution,  though  less  significant  than  in  the  negative  electrode 
for  two  reasons.  First,  at  the  onset  of  the  30  A  discharge, 


Fig.  8.  Distribution  of  active  material  surface  concentration  across  cell  for  first 
30  s  of  HPPC  test  at  58.3%  SOC. 

initial  reaction  distribution  is  already  more  uniform  in  the 
positive  electrode.  Second,  and  more  significant,  the  positive 
electrode  open-circuit  potential  function  is  almost  seven  times 
more  sensitive  to  changes  in  concentration  than  the  negative 
electrode  function  at  58.3%  SOC.  Small  changes  in  positive 
electrode  active  material  surface  concentration  significantly 
penalize  reaction,  and  for  this  reason,  redistribution  of  reaction 
due  to  solid  diffusion  limitations  occurs  much  quicker  in  that 
electrode. 

Redistribution  of  Li  also  occurs  upon  cell  relaxation  at  the 
end  of  the  18  s-long  30  A  discharge.  Fig.  7  shows  a  second  step 
change  in  transfer  current  density  around  19  s,  appearing  qual¬ 
itatively  as  the  mirror  image  of  the  step  change  at  1  s  when 
the  galvanostatic  load  was  first  applied.  At  locations  near  the 
separator,  a  recharging  process  begins  while  at  locations  near 
the  current  collector,  discharge  continues  even  after  the  load  is 
removed.  The  process  is  driven  by  the  gradient  in  local  equilib¬ 
rium  potential,  dU/dx  (directly  related  to  the  solid  phase  surface 
concentration  gradient,  dcs^/dx)  and  continues  until  surface  con¬ 
centrations,  cs?e,  are  once  again  evenly  distributed.  During  this 
redistribution  process  the  net  balance  of  reaction  across  each 
electrode  is  zero. 

Relaxation  reaction  redistribution  is  less  significant  in  the 
positive  electrode,  where  only  a  minimal  solid  phase  surface 
concentration  gradient,  3cs?e/ dx,  arose  during  the  30  A  discharge. 
The  process  lasts  on  the  order  of  10  s,  compared  to  several  min¬ 
utes  for  the  negative  electrode.  In  both  electrodes,  solid  phase 
bulk  concentrations  rise  and  fall  at  roughly  the  same  rate  as 
surface  concentrations  throughout  the  redistribution  process, 
indicating  that  the  time  scale  of  reaction  redistribution  is  much 
faster  than  solid  phase  diffusion.  Localized  concentration  gra¬ 
dients  within  individual  solid  particles  (from  bulk  to  surface), 
3 cs/3r,  relax  so  slowly  that  the  63  s  long  pulse  power  test  amounts 
to  little  more  than  a  transient  discharge  and  charge  on  the  sur¬ 
face  of  the  active  material  particles  with  inner  bulk  regions 
unaffected. 

4.2.  Rate  capability 

Fig.  9  presents  model-predicted  constant  current  discharge 
capability  from  50%  SOC.  Shown  in  the  bottom  window  of 
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Fig.  9.  Solid  phase  surface  concentration  (top),  minimum  electrolyte  concen¬ 
tration  (middle),  and  time  (bottom)  at  end  of  galvanostatic  discharge  for  rates 
from  10  to  40  C. 


Fig.  10.  Dimensionless  distribution  of  concentration  within  an  active  material 
particle  at  various  times  during  galvanostatic  (dis)charge. 


Fig.  9,  a  40  C  rate  current  (240  A)  can  be  sustained  for  just 
over  6  s  before  voltage  decays  to  the  2.7  V  minimum.  The  top 
window  of  Fig.  9  shows  active  material  surface  concentration 
in  the  negative  (left  axis)  and  positive  (right  axis)  electrodes 
at  the  end  of  discharge  across  the  range  of  discharge  rates. 
Electrode-averaged  rather  than  local  values  of  surface  concen¬ 
tration  are  presented  to  simplify  the  discussion.  The  electrodes 
are  fairly  well  balanced,  indicated  by  end  of  discharge  surface 
concentrations  near  depletion  and  saturation  in  the  negative  and 
positive  electrodes,  respectively.  End  of  discharge  voltage  is 
predominantly  negative  electrode-limited  as  stoichiometries  of 
v  =  cS5e/cs,max-  <0.05  causes  a  rapid  rise  in  U-.  Surface  active 
material  utilization  decreases  slightly  with  increasing  C-rate  due 
to  increased  ohmic  voltage  drop. 

In  the  present  model,  electrolyte  Li+  transport  is  sufficiently 
fast  that  electrolyte  depletion  does  not  play  a  limiting  role  at 
any  discharge  rate  from  50%  SOC.  In  the  worst  case  of  30  C, 
the  minimum  value  of  local  electrolyte  concentration,  occurring 
at  the  positive  electrode/current  collector  interface,  is  around 
50%  of  average  concentration,  ce,o-  For  rates  less  than  30  C, 
the  reduced  current  level  results  in  lesser  electrolyte  concen¬ 
tration  gradients,  while  for  rates  greater  than  30  C,  the  shorter 
duration  of  discharge  time  results  in  a  smaller  concentration 
gradient  at  end  of  discharge.  If  we  induce  sluggish  diffusion 
by  reducing  De  (a  similar  effect  may  be  induced  in  cell  design 
by  reducing  porosity),  electrolyte  concentration  in  the  positive 
electrode  comes  closer  to  depletion  with  the  worst  case  mini¬ 
mum  value  of  ce  occurring  at  lesser  current  rates.  Lowering  De 
by  one  order  of  magnitude  for  example,  results  in  a  battery  lim¬ 
ited  in  the  10-20  C  range  by  electrolyte  phase  transport,  with 
higher  and  lower  current  rates  still  controlled  by  solid  phase 
transport. 


4.3.  Solid-state  diffusion  limited  current 

Under  solid  phase  transport  limitations,  simple  relationships 
may  be  derived  to  predict  maximum  current  available  for  a  given 


pulse  time.  Substituting  dimensionless  variables: 


r  = 


R, 


x  = 


Dst 
R2  ’ 


cs(r,  r)  = 


cs(r,  r)  -  cs,o 


c 


s,max 


r  = 


ju  /?< 


DsasFcsmax 


(11) 


into  Eq.  (4)  yields  the  dimensionless  governing  equation: 


dcs  _  1  d  f  2dcs 
dr  r2  dr  \  dr 


(12) 


with  initial  condition  cs(r,  r  =  0)  =  OVr  and  boundary  condi¬ 
tions: 


dcs 

dr 


=  0 , 


r= 0 


dCs 

dr 


=  JU 


(13) 


r=l 


The  solution  given  by  Carslaw  and  Jaeger  [25]  is 


cs(r,  r)=  j 


-•Li 


1  9 

3rH - (5r2 

10 


oo 


_  ^_2y^sin(A„r)  exp(-^r) 


n  —  1 


X2n  sin(A„) 


(14) 


where  the  eigenvalues  are  roots  of  kn  =  tan(A.n).  Fig.  10  shows 
distribution  of  Li  concentration  along  the  radius  of  an  active 
material  particle  during  galvanostatic  discharge  or  charge  for 
dimensionless  times  ranging  from  r=  10-6  to  10_1  which,  for 
reference,  correspond  to  current  pulses  lasting  0.05-500  s  using 
negative  electrode  parameters  from  Table  2. 

Surface  concentration  and  depth  of  penetration  into  the  active 
material,  both  of  practical  interest  for  HEV  pulse-type  operation, 
can  be  obtained  from  Eq.  (14).  Surface  concentration,  shown  for 
the  present  model  to  cause  end  of  discharge  as  the  negative  elec¬ 
trode  nears  depletion,  is  calculated  by  evaluating  Eq.  (14)  at 
r  =  1 .  Penetration  depth,  8 ,  providing  a  measure  of  active  mate¬ 
rial  accessible  for  short  duration  pulse  events,  is  calculated  by 
finding  the  point  along  the  radius  where  the  concentration  profile 
is  more  or  less  equal  to  the  initial  condition.  A  99%  penetration 
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Table  3 

Empirical  formulae  fit  to  solid-state  diffusion  PDE  exact  solution,  Eq.  (14) 


1  %  error  bounds 


Dimensionless  surface  concentration 

j!’6  =  1.139  Vr 

ju 

(16) 

0  <  T  <  1  x  10“4 

lls;e  =  l.myi  1.25r 
ju 

(17) 

0<t<8  x  10~2 

Dimensionless  99%  penetration  depth 

8  =  3.24/r 

(18) 

0  <  r  <  1  x  10-3 

8  =  3.23Vr+  1-89t 

(19) 

0<t<2  x  10-2 

depth,  8  =  Rs  —  r,  is  defined  using  the  location  r  resulting  in  a 
root  of  the  formula: 


Cs(r,  r) 

cs(l,  O 


(15) 


with  a  =  0.99.  Expressed  as  a  fraction  of  total  radius,  dimension¬ 
less  penetration  depth,  8  =  8/ Rs,  is  a  function  of  dimensionless 
time  only. 

Empirical  expressions  for  dimensionless  surface  concentra¬ 
tion,  cs>e,  and  dimensionless  penetration  depth,  5,  are  fit  to  the 
results  of  Eq.  (14)  and  presented  in  Table  3.  Functions  of  the 
form  f{x)  =  Ca/t  provide  good  resolution  at  short  times  of 
x  <  10~3,  corresponding  in  our  model  to  pulses  lasting  fewer 
than  ~5  s.  Resolution  may  be  extended  one  to  two  orders  of  mag¬ 
nitude  in  x  by  using  functions  of  the  form  fix)  =  C^fx  +  Dx. 
The  latter  functions  are  plotted  versus  the  exact  solution  (14)  in 
Fig.  11. 

Eq.  (17)  in  Table  3  may  be  used  in  lieu  of  the  present  elec¬ 
trochemical  model  to  predict  negative  electrode  surface  concen¬ 
trations  for  pulses  shorter  than  400  s,  or,  conversely,  to  predict 
limiting  currents  at  rates  greater  than  ^9  C  caused  by  depleted 
active  material  surface  concentration  at  end  of  discharge.  By 
combining  Eqs.  (11)  and  (17)  under  the  assumption  of  uniform 
current  density,  =  I/AL-,  we  obtain  an  empirical  relation- 


X 


Fig.  11.  Dimensionless  active  material  surface  concentration  and  penetration 
depth  versus  time  (top).  Percent  error  in  empirical  correlations  of  the  form 
f(x)  =  C-s/t  +  Dr  (bottom). 


ship  for  surface  concentration  as  a  function  of  current  and  time: 
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(20) 


valid  for  t  <  0.08 Ds/R^.  Alternatively,  given  initial  stoichiom¬ 
etry,  vo,  and  surface  stoichiometry  at  end  of  discharge,  vs,e  final? 
the  maximum  current  available  for  a  pulse  discharge  lasting  t 
seconds  will  be 


Anax  —  (tq  3Cs?e  final) 


L-Aa*-Fc 


s,max- 


1.122 


\ft  + 


1.25 

R*- 


(21) 


While  the  theoretical  maximum  current  will  be  obtained  under 
the  condition  vs?e  final  =  0,  i.e.  complete  surface  depletion,  Fig.  9 
showed  the  present  model  to  exhibit  an  end  of  discharge  surface 
stoichiometry  around  0.03,  with  some  rate  dependency.  Under 
uniform  initial  conditions,  the  initial  stoichiometry,  xq,  is  sim¬ 
ply  a  function  of  SOC.  For  a  recently  charged  or  discharged 
battery  with  nonuniform  initial  concentration,  a  better  predic¬ 
tion  of  maximum  pulse  current  may  be  obtained  by  replacing  vo 
in  Eq.  (21)  with  a  stoichiometry  averaged  across  the  penetration 
depth  or  “pulse-accessible”  region  [26]. 


4.4.  Suitability  of  solid-state  diffusion  approximations  for 
cell  modeling 


Introduced  in  Section  2  and  detailed  in  Appendix  A,  the 
present  work  utilizes  a  fifth-order  finite  element  approximation 
for  solid-state  diffusion  (4)  and  incorporates  that  submodel  into 
the  ID  electrochemical  model  as  a  finite  difference  equation, 
that  is,  local  values  of  cs?e  are  calculated  using  values  of  cs?e  and 
jLl  from  the  previous  five  time  steps: 


.[*]  =  f(\rVk-n  r[k-5h  v;U[k] 
s  e  J  fLcs,e  ?  •  •  •  ?  cs,e  J’  L J  ? 


•  • ,  7Li[*-5]]) 


(22) 


Here  we  compare  the  finite  element  submodel  to  the  analyti¬ 
cal  solution  employed  by  Doyle  et  al.  [11]  and  the  polynomial 
profile  model  of  Wang  et  al.  [19].  Comparisons  are  made  in 
the  frequency  domain  in  order  to  remove  the  influence  of  a 
particular  type  of  input  (pulse  current,  current  step,  constant  cur¬ 
rent,  etc.)  by  taking  the  Laplace  transform  of  each  time  domain 
model,  expressing  the  input/output  relationship  as  a  transfer 
function  in  the  Laplace  variable  s,  and  substituting  s=jco  to 
calculate  the  complex  impedance  at  frequency  co.  A  capital¬ 
ized  variable  denotes  that  variable’s  Laplace  transform,  that 
is,  CSje(^)  =  E{cs,e(0}?  and  an  overbar  denotes  a  dimensionless 
variable.  Define 


CsAs)  = 


Q,eW  cs,0 


C 


s,max 


CO  =  CO 


D, 


jLi(s )  =  yLi(s) 


OsFDtCs, 


max 


(23) 


In  the  Laplace  domain,  a  compact  analytical  solution  to  Eq. 
(4)  is  readily  available.  The  exact  transfer  function  expressing 
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dimensionless  surface  concentration  versus  dimensionless  reac¬ 
tion  current  given  by  Jacobsen  and  West  [28]  is 


QeO)  _  tanh(^) 
7Ll0)  tanh(^)  —  x//  ’ 


(24) 


where  ^  =  Rs\/s/Ds. 

Doyle  et  al.  [11]  provide  two  analytical  series  solutions  in 
the  time  domain,  one  for  short  times  and  one  for  long  times.  In 
Appendix  B,  we  manipulate  Doyle’s  formulae  to  arrive  at  the 
short  time  transfer  function: 


CsAs) 

JLi(s) 


oo 


i  -1 


1  —  \[r  +  2V^y^exp(— 2ml/) 


n  —  1 


(25) 


and  the  long  time  transfer  function: 


CsAs) 

Ju(s ) 


(26) 


The  frequency  response  (magnitude  and  phase  angle)  of  trun¬ 
cated  versions  of  the  short  and  long  time  transfer  functions  are 
compared  to  the  exact  transfer  function  (24)  in  Fig.  12,  showing 
the  short  time  solution  to  provide  good  agreement  at  high  fre¬ 
quencies  and  the  long  time  solution  at  low  frequencies.  Note  that 
the  short  time  transfer  function  does  not  change  much  beyond  the 
first  term  of  the  series.  A  good  strategy  to  piece  together  Doyle’s 
two  solutions  is  to  use  one  term  of  the  short  time  solution  for 
r  =  Dst/R2s  <  0.1  (corresponding  to  co  =  6  x  101  in  Fig.  12) 
and  around  100  terms  of  the  long  time  solution  for  r  >  0.1. 

Reaction  current  appears  in  Eq.  (4)  as  a  time  depen¬ 
dent  boundary  condition  which  Doyle  accommodates  using 
a  Duhamel  superposition  integral.  Numerical  solution  of  this 
convolution-type  integral  requires  that  a  time  history  of  all  pre¬ 
vious  step  changes  in  surface  concentration  be  held  in  memory 
and  called  upon  at  each  time  step  to  reevaluate  the  integral.  So 
while  the  analytical  solution  is  inarguably  the  most  accurate 
approach,  it  can  be  expensive  in  terms  of  memory  and  computa¬ 
tional  requirements,  particularly  in  situations  requiring  a  small 


Fig.  12.  Frequency  response  of  short  and  long  time  analytical  solutions  used 
in  Ref.  [11]  for  solid-state  diffusion  in  spherical  particles,  compared  to  exact 
frequency  response  from  Ref.  [28]. 


Fig.  13.  Frequency  response  of  parabolic  profile  solid-state  diffusion  submodel 
from  Ref.  [19]  and  fifth-order  finite  element  solid-state  diffusion  submodel  (used 
in  this  work),  compared  to  exact  frequency  response  from  Ref.  [28]. 


time  step  but  long  simulation  time  (driving  cycle  simulations,  for 
instance)  or  in  situations  requiring  a  large  grid  mesh  (2D  or  3D 
simulations  incorporating  realistic  cell  geometry,  for  instance). 

Approximate  solutions  to  Eq.  (4)  are  appropriate  so  long  as 
they  capture  solid-state  diffusion  dynamics  sufficiently  fast  for 
a  particular  investigation.  Wang  et  al.  [19]  assume  the  concen¬ 
tration  profile  within  the  spherical  particle  is  described  by  a 
parabolic  profile  cs(r,  t)=A(t)  +  B(t)r2,  and  thus  formulate  a 
solid-state  diffusion  submodel  which  correctly  captures  bulk 
dynamics  and  steady  state  concentration  gradient,  but  otherwise 
neglects  diffusion  dynamics.  Derived  in  Appendix  C,  the  trans¬ 
fer  function  of  the  parabolic  profile,  or  steady  state  diffusion, 
model  is 


CsAs)  3  1 

JLl(s)  \j/2  5 


(27) 


Shown  in  Fig.  13  versus  the  exact  transfer  function  (24),  the 
parabolic  profile  model  is  valid  for  low  frequencies,  6b  <  10, 
or  long  times,  r  =  Dst/R 2  >  0.6.  Substituting  values  from  the 
present  model’s  negative  electrode  (Ds_  =  2.0  x  10  12  cm2  s  1 , 
Rs_  =  1.0  x  10-4  cm),  the  parabolic  profile  model  would  cor¬ 
rectly  predict  surface  concentration  only  at  times  longer  than 
3000  s.  For  electrochemical  cells  with  sluggish  solid-state  diffu¬ 
sion,  the  parabolic  profile  model  will  correctly  capture  low-rate 
end  of  discharge  behavior,  but  is  generally  inappropriate  in  the 
modeling  of  high  rate  (>2  C)  or  pulse  type  applications  [27] . 

We  find  spatial  discretization  of  Eq.  (4)  yields  low-order 
solid-state  diffusion  models  with  more  accurate  short  time  pre¬ 
diction  compared  to  polynomial  profile  models  [20].  Recasting 
the  fifth- order  finite  element  model  from  Appendix  A  in  nondi- 
mensional  form: 


CsAs) 

Ju(s) 


b\S5  +  +  £>3  S3  +  b^S2  +  b$s  + 

a\s5  +  a2S4  +  a^s2  +  <24  s2  +  a$s  +  a 6 

(28) 


Fig.  13  shows  the  present  model  to  provide  good  approximation 
of  the  exact  transfer  function  (26)  for  60  <  105 ,  and  thus  be  valid 
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for  dimensionless  times  r  >  6  x  10~5  (or  t>  0.3  s  for  the  present 
model’s  negative  electrode).  Regardless  of  what  solution  tech¬ 
nique  is  employed  for  solid-state  diffusion  in  an  electrochemical 
cell  model,  if  the  objective  is  to  match  high-rate  (~40  C)  pulse 
behavior  and  predict  transport  limitations  on  a  short  (~5  s)  time 
scale,  that  technique  must  be  valid  at  very  short  times. 

5.  Conclusions 


sented  as  ODEs  in  state  space  form: 


G,  1 

"G,r 

cs,l 

Os, 2 

=  A 

cs,2 

+  ByLl. 

G,e  ~  C 

cs,2 

+  i)/Ll. 

-6s,n  _ 

_  cs,n  _ 

_  Os,n  _ 

A  fifth-order  finite  element  model  for  transient  solid-state 
diffusion  is  incorporated  into  a  previously  developed  ID  electro¬ 
chemical  model  and  used  to  describe  low-rate  constant  current, 
hybrid  pulse  power  characterization,  and  transient  driving  cycle 
data  sets  from  a  lithium  ion  HEV  battery.  HEV  battery  models  in 
particular  must  accurately  resolve  active  material  surface  con¬ 
centration  at  very  short  dimensionless  times.  Requirements  for 
the  present  model  are  t  =  Dst/R 3  ~  10-3  to  predict  40  C  rate 
capability  and  r~2x  10~5  to  match  current/voltage  dynamics 
at  10  Hz. 

Dependent  on  cell  design  and  operating  condition,  end  of 
pulse  discharge  may  be  caused  by  negative  electrode  solid  phase 
Li  depletion,  positive  electrode  solid  phase  Li  saturation,  or 
electrolyte  phase  Li  depletion.  Simple  expressions  developed 
here  for  solid-state  diffusion-limited  current,  applicable  in  either 
electrode,  may  aid  in  the  interpretation  of  high-rate  experimental 
data.  While  the  present  work  helps  to  extend  existing  litera¬ 
ture  into  the  dynamic  operating  regime  of  HEV  batteries,  future 
work  remains  to  fully  characterize  an  HEV  battery  in  the  lab¬ 
oratory  and  develop  a  fundamental  model  capable  of  matching 
current/voltage  data  at  very  high  rates. 


where  the  n  states  of  the  system  are  the  radially  distributed  val¬ 
ues  of  concentration  cs,i, . . .,  cs?n,  at  finite  element  node  points 
1 For  the  linear  PDE  (4)  with  constant  diffusion  coeffi¬ 
cient,  Ds,  the  matrix  A  is  constant  and  tri-diagonal. 

The  linear  state  space  system  (A.l)  can  also  be  expressed  as 
a  transfer  function: 

cs  Q(s)  b\sn  +  b2Sn~l  +  •  •  •  +  bn-\s  +  bn 

%  G(s)  =  — — — - f- -  n  1  n  (A. 2) 

jLl(,s)  a\sn  +  ci2Sn  1  +  •  •  •  +  an-\s  +  an 

with  constant  coefficients  at  and  bi  [30].  While  either  (A.l)  or 
(A. 2)  could  be  numerically  implemented  using  an  iterative  solu¬ 
tion  method,  we  exploit  the  linear  structure  of  (4)  and  express  the 
system  as  a  finite  difference  equation  with  explicit  solution.  To 
discretize  (A. 2)  with  respect  to  time,  we  perform  a  z-transform 
using  Tustin’s  method: 

GT(z)  =  G(s)\s=(2/ts)((z-1)/(z+1))’  (A'3) 

resulting  in  an  nth-order  discrete  transfer  function: 

Cs,e(z)  r-  t  \  h\Zn  +  h2Zn~l  H - \-hn-iz  +  hn 

jU(z)  gizn  +  g2Zn  1  H - 1-  gn-iz  +  gn 

(A.4) 
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Appendix  A 

A.L  Solid-state  diffusion  finite  element  model 

The  transient  phenomenon  of  solid-state  Li  diffusion  is  incor¬ 
porated  into  the  previously  developed  macroscopic  model  of  Gu 
and  Wang  [3].  While  the  governing  Eq.  (4)  describes  solid  phase 
concentration  along  the  radius  of  each  spherical  particle  of  active 
material,  the  macroscopic  model  requires  only  the  concentration 
at  the  surface,  cs,e(0>  as  a  function  of  the  time  history  of  local 
reaction  current  density,  jL\t). 

We  transform  the  PDE,  Eq.  (4),  from  spherical  to  planar 
coordinates  using  the  substitution  v(r)  =  rcs(r)  [28,29]  and  dis¬ 
cretize  the  transformed  equation  in  the  r-direction  with  n  linear 
elements.  (The  present  model  uses  five  elements  with  node 
points  placed  at  {0.7,0.91,0.97,0.99,1.0}  x  Rs.)  Transformed 
back  to  spherical  coordinates,  the  discretized  system  is  repre- 


with  constant  coefficients  hi  and  gi.  Computation  is  thus  reduced 
to  an  explicit  algebraic  formula  with  minimal  memory  require¬ 
ments.  Solution  for  cs?e  requires  that  local  values  of  cs?e  and  jLl 
be  held  from  only  the  previous  n—  1  time  steps. 


Appendix  B 


B.l.  Transfer  functions  of  short  and  long  time  solid-state 
diffusion  analytical  solutions 


Doyle  et  al.  [11]  employ  an  analytical  solution  to  Eq.  (4)  and 
embed  it  inside  a  Duhamel  superposition  integral  to  accommo¬ 
date  the  time  dependent  boundary  condition.  They  provide  two 
integral  expressions  for  the  response  of  reaction  current,  jLl(r), 
to  a  step  in  surface  concentration,  A cs,e,  at  r  =  0,  each  in  the 
form 


1  Rs  1 

asF  Ds  Acs?e 


jU(0  d£ 


(B.l) 


The  short  time  expression  (Eq.  (B. 6)  of  [11])  is 


(B.2) 
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while  the  long  time  expression  (Eq.  (B.5)  of  [1 1])  is 


2  °°  1  99 
a(  r)  =  -^-[1  exp(-^27r2r)]. 

-rrA  /  ■</  mZ 


(B.3) 


7i  — '  n 

n= 1 


Little  detail  is  given  in  the  derivation  of  these  expressions.  Dif¬ 
ferentiating  Eqs.  (B.1)-(B.3)  with  respect  to  r  and  solving  for 
jLi(t),  we  recover  short  time  solution: 


t  i  Ds 

jL  (r)  =asF— 

Ac 


oo 


1  - 


1  +  ( - 


rv 


JVC 


n  —  1 


Ac 


s,e 

(B.4) 


and 

cs,eO  -  cs,avgw  =  -A—  yLio,  (c.4) 

which  when  combined  to  eliminate  Cs?avg,  provides  the  transfer 
function: 

_  1 

JLl(s)  asF  Ds 

Substituting  dimensionless  variables  JL\  CSiQ,  and  x/z  into  Eq. 
(C.5)  yields  Eq.  (27),  used  in  Section  4.4. 

References 


3  Ds  1 

2M  +  5 


(C.5) 


and  long  time  solution: 


r  ;  Ds 

jL  (r)  =asF— 

As 


OO 


— 2y^exp(— n2n2r) 


n  —  1 


Ac 


s,e 


(B.5) 


no  longer  in  integral  form.  Taking  the  Laplace  transform  of  Eqs. 
(B.4)  and  (B.5)  and  recognizing  that  the  transform  of  the  step 
input  is  Cs,e(s)  =  Acs,e/s,  we  find  the  short  time  transfer  function: 


JU(s)  Ds 

=  asF 


C S,e(s) 


R < 


1  -  R< 


D 


+  2  Rf 


xexp  [  —2nRs 


D 


and  the  long  time  transfer  function: 


(B.6) 


JU(s)  D s 

=  asF 


C s,eC) 


R< 


OO 


-2E 


«=1 


s  +  n2  7t2 


(B.7) 


Taking  the  reciprocal  of  Eqs.  (B.6)  and  (B.7)  and  substituting 
dimensionless  variables  7Ll,  and  x/z  yields  Eqs.  (25)  and 
(26)  respectively,  used  in  Section  4.4. 


Appendix  C 

C.l.  Transfer  function  of  parabolic  profile  solid-state 
diffusion  approximate  model 


Wang  et  al.  [19]  assume  concentration  distribution  within  a 
spherical  active  material  particle  to  be  described  by  a  parabolic 
profile.  Integrating  the  two  parameter  polynomial  with  respect 
to  Eq.  (4),  they  reduce  the  problem  of  determining  surface  con¬ 
centration,  cs?e(0>  as  a  function  of  reaction  current,  jLl(t),  down 
to  the  solution  of  one  ODE: 


9cs,avg 

dt 


3 

asFRs 


and  one  interfacial  balance: 


(C.l) 


cs,e  cs,avg 


R < 


5  a*FD, 


J 


:Li 


(C.2) 


Taking  the  Laplace  transform  of  Eqs.  (C.l)  and  (C.2)  yields: 


^C^aygCs’)  — 


a*FR, 


Ju{s) 


(C.3) 
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